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1. Introduction 

Dark matter (DM) distributes in the Universe in halos that host galaxy clusters, galaxies and 
galactic DM “clumps”. A promising way to identify the nature of DM and measure its properties is 
to search for the Standard Model (SM) particles produced in its annihilation or decay at these sites. 
Gamma rays and neutrinos are ideal messengers for directional DM searches, since they are the 
only stable neutral SM particles, and can thus travel from their production sites to Earth unaffected 
by magnetic deflection. 

Current gamma-ray instruments like the Fermi-LP^ in space, the ground-based Cherenkov 
telescopes MAGIC, VERITAS and H.E.S.S., and the new-generation water Cherenkov detector 
HAWC, as well as neutrino telescopes like IceCube and Antares, are sensitive to overlapping and 
complementary DM particle mass ranges (from ~1 GeV to ~100 TeV). All of these instruments 
have dedicated programs to look for DM signals coming from, e.g., the Galactic center and halo 
[1, 2, 3, 4, 5, 6, 7] , galaxy clusters [8, 9, 10, 11, 12], or dwarf spheroidal satellite galaxies (dSphs) 
of the Milky Way (MW) [13, 14, 15, 16]. 

The universality of DM properties allows the combination of data from different experiments 
and/or observational targets into a global, sensitivity-optimized search [17]. Eor a given DM par¬ 
ticle model, a joint likelihood function can be written as the product of the partial contributions 
from each of the measurements/instruments. The advantage of such an approach is that the details 
of each experiment do not need to be combined or averaged. We have implemented this analysis 
framework, applicable to observations from gamma-ray and neutrino instruments, and applied it to 
the MAGIC and Fermi-EAT observations of dSphs. 

The MW dSphs are associated to the Galactic DM sub-haloes, predicted by N-body cosmolog¬ 
ical simulations, that have attracted enough baryonic mass to start stellar activity (other sub-halos 
may remain completely dark). MW dSphs have very high mass-to-light ratios, being the most 
DM-dominated systems known so far [18]. MW dSphs also have the advantage of being free of 
astrophysical gamma-ray sources and, while they are relatively close, they still appear as quasi- 
point-like sources for gamma-ray and neutrino telescopes, with relatively high expected fluxes. In 
addition, the possibility of determining the DM distribution and its uncertainty using a common 
methodology [19] allows a straightforward combination of their observations into a sensitivity- 
optimized global analysis. 

In this paper, we present our new global analysis framework and the results of applying it to 
MAGIC and F?rnj/-EAT observations of dSphs. 

2. MAGIC and Fermz-LAT data samples 

The Florian Goebel MAGIC telescopes are located at the Roque de los Muchachos Observa¬ 
tory (28.8° N, 17.9° W; 2200 m above sea level), at the Canary Island of Ea Palma (Spain). MAGIC 
is a system of two telescopes that detect Cherenkov light produced by the atmospheric showers ini¬ 
tiated by cosmic particles entering the Earth atmosphere. Cherenkov images of the showers are 
projected by MAGIC reflectors onto the photo-multiplier tube (PMT) cameras, and are used to 
reconstruct the calorimetric and spatial properties of the primary particle, as well as for its iden¬ 
tification. Thanks to its large reflectors (17 meter diameter), plus its high-quantum-efficiency and 
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low-noise PMTs, MAGIC achieves high sensitivity to Cherenkov light, hence low energy thresh¬ 
old. The MAGIC telescopes are able to detect cosmic gamma rays in the very-high-energy domain, 
i.e. in the range between ~50 GeV and ~50 TeV. 

For our study, we use MAGIC data corresponding to 158 hours of observations of Segue 1 [13], 
the deepest observations of any dSph by any Cherenkov telescope. The data were taken between 
2011 and 2013 [20]. 

The Fermi-LKT is a pair-conversion telescope sensitive to gamma rays in the range from 
20 MeV to more than 300 GeV [21]. With its large field of view (2.4 sr), the Fenni-LAI is able 
to efficiently survey the entire sky. Indeed, since its launch in August 2008, the Fermi-LAI has 
primarily operated in a survey mode that scans the entire sky every 3 hours. The survey-mode 
exposure coverage is fairly uniform over the sky with variations of at most 30% with respect to the 
average exposure. The Femii-LPJ' source sensitivity, which is limited by the intensity of diffuse 
backgrounds, shows larger variations but is relatively constant at high galactic latitudes (b > 10°). 

In this work, we use a Fermi-LAI data sample corresponding to 6 years of observations of 
15 dSphs, processed with the latest (Pass 8) data analysis [16]. Events are selected with energies 
between 500 MeV and 500 GeV in a 10° x 10° region of interest (ROI) centered on each dSph. 
EAT likelihoods for a given DM model are constructed from the bin-by-bin likelihoods of [16] and 
do not involve any reanalysis of the EAT photon data. 


3. Analysis 


3.1 Dark Matter annihilation flux 


The gamma-ray (or neutrino) flux produced by DM annihilation in a given largel and defecfable 
af Earth by an insfrumenf observing a field of view AQ is given by: 


dE ' 471 dE 


(3.1) 


where widm is the mass of the DM particle, (av) the thermally-averaged annihilation cross section, 
dN/dE the average gamma-ray spectrum per annihilation reaction (for neutrino this term includes 
the oscillation probability between target and Earth), and 

7(AD)= [ dH [ dlp^{l,a) (3.2) 

JaQ. Jl.o.s. 

is the so-called astrophysical factor (or simply J-factor), with p being the DM density, and the 
integrals running over AD and the line of sight (l.o.s.), respectively. 

Using PYTHIA simulation package version 8.205 [22], we have computed the average gamma- 
ray spectra per annihilation process {dN/dE) for a set of DM particles of masses between 10 GeV 
and 100 TeV, annihilating into SM pairs bb and T+T^. Eor each channel and mass, we average the 
gamma-ray spectrum resulting from 10^ decay events of a generic resonance with mass 2 x /mdm 
into the specified pairs. Eor each simulafed evenf, we frace all fhe decay chains, including fhe muon 
radiative decay (p^ —)• nof active in PYTHIA by defaull), down fo sfable particles. 

Eor fhe compufafion of fhe J-faclors we follow fhe approach by Martinez [19]. The disfribufion 
of DM in fhe halos of fhe differenl dSphs are parameferized following a Navarro-Erenk-While 
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profile (NFW) [23]: 


p(r) = 


Pori 


(3.3) 


r(r, + r)2 

where and po are the NFW scale radius and characteristic density, respectively, and are deter¬ 
mined from a fit to the dSph stellar density and velocity dispersion profiles. The properties of the 
dSphs used in our analysis, including the J-factors and their uncertainties, can be found in [16]. 

We use templates for the DM emission in each dSph normalized to its J-factor integrated to a 
radius of 0.5° from the halo center (Tobs.i)- The 0.5° integration region encompasses more than 90% 
of the annihilation flux for our dSph halo models which have halo scale radii between 0.1° and 0.4°. 
In the LAT analysis the intensity templates are used to construct a three-dimensional model for the 
expected DM signal as a function of space and energy within the 10° x 10° ROI centered on each 
dSph. The MAGIC analysis uses a one-dimesional likelihood for the photon energy distribution 
within a signal aperture of radius 0.122°. The observable flux in the MAGIC analysis of Segue 1 is 
therefore reduced by a factor of ~ 1.6 with respect to the LAT analysis of the same target. However 
we note that because the signal aperture is matched to the angular size of the MAGIC PSF, this 
truncation has a negligible impact on the sensitivity of the analysis. 


3.2 Likelihood analysis 


For each considered annihilation channel and DM particle mass, we compute the profile like¬ 
lihood ratio as a function of (cJv): 


Ap((av)|^) 




(3.4) 


with & representing the data samples and V the nuisance parameters, (ov) and V are the values 
maximizing the joint likelihood function {^), and V the value that maximizes ^ for a given value 
of {ov). The likelihood function can be written as: 


^((av);v|^) 


Mai'get 

n 


1=1 


.^i{{ovy,Ji,pii I &i) ■ /{Ji 17obs,0 Oi) , 


(3.5) 


with the index i running over the different targets (dSphs in our case); 7, is the J-factor for the 
corresponding target (see Equation 3.2); ju, denotes any additional nuisance parameters; and the 
target-related input data. ^ is the likelihood for T,-, given measured logjQ(7obs,i) and its uncertainty 
Oi [16], i.e.: 


^ (7; I 7obs,i) Oi ) 


ln( 10)7obs,i 


(3.6) 


The likelihood function for a particular target (.it]) can in turn be written as the product of the 
likelihoods for different instruments (represented by the index j), i.e.: 


^in strum ent 

.^i{{ovy,ji,Hi I &i) = n .^ij{{ov)-Ji,Hij I %) 

y=i 


(3.7) 


where and represent the nuisance parameters and input data sample for the given target i 
and instrument j. 


4 






MAGIC/Fermi-LAT combined dark matter searches 


Matthew Wood 


Equations 3.4, 3.5 and 3.7 are generic, i.e. they are valid for any set of instruments and ob¬ 
served targets. In addition, they allow merging of the results from different instruments and targets, 
starting from tabulated values of vs. (av) for a fixed value of 7, and profiled with respect to 
ju,y. These values can be produced and shared by the different experiments without the need of 
releasing or sharing any of the internal information used to produce them. 

For this work, we use the Fenni-LKV likelihood values vs. energy flux, tabulated for different 
targets and in energy bins, and released by the Fenni-LKA Collaboration [16]. In the case of 
MAGIC, the likelihood is obtained following the method described in Refs. [17] and [13]. 

4. Results 

We compute one-sided, 95% confidence level upper limits to (av) by numerically solving the 
equation —21nAp((av)2,7i | &) = 2.71, for {ov)2.i\- We consider Xp restricted to non-negative 
(av) values only. 

Figure 1 shows the 95% confidence level limifs to (av) for DM particles with masses between 
10 GeV to 100 TeV annihilating into SM pairs {bb and t+t^) obtained from the combination 
of Fermi-LAJi (15 dSphs and Segue 1 alone) and MAGIC Segue 1 observations. The 68% and 
95% containment bands are computed from the distribution of upper limits obtained from 300 
analysis realizations consisting of Fermi-LA^ observations of empty fields of view combined wifh 
MAGIC fasf simulations (assuming in both cases equal exposures as for the real data), and J-factors 
randomly selected according to the probability density function (PDF) in Equation 3.6. The blank 
fields constifufing fhe FAT realizations were selected by choosing randomy sky positions with 
\b\ > 30° centered at least 0.5° from a source in the 3FGF catalog. MAGIC fast simulations consist 
of a set of event energies randomly generated from the background PDF (see Ref. [13] for details) 
for both signal and background regions. 

We find no positive signal of DM in our data sample. As expected, limits in the low and high 
ends of the considered mass range are dominated by Fermi-LA^ and MAGIC observations, respec¬ 
tively, where the combined limits coincide with the individual ones. The combination provides a 
significant improvement in the range between ~1 and ~30 TeV (for bb) or ~0.2 and ~2 TeV (for 
T+T^), with a maximum improvement of the combined limits with respect to the individual ones 
by a factor ~2 at a mass of 500 GeV (for bb) and 3 TeV (for t+t^). 

MAGIC individual results shown here may differ from those presented in Ref. [13] by up 
to a factor ~4, which needs a dedicated explanation. First, we note that the data, instrument re¬ 
sponse functions, and likelihood functions are identical in both works. Aside from enlarging the 
explored DM mass range, and in order to homogenize MAGIC and Fermi-FA^ analyses, we have 
introduced the following differences between the two works: i) the J-factor for MAGIC cut of 
7(0.122°) = 2.2 X 10^^ GeV^ cm^^(following Ref. [19] and assuming an NFW DM density pro¬ 
file); ii) include the statistical uncertainties in the determination of the J-factor; Hi) use the Fermi- 
FAT prescription for limits close to bounds of the physical region ((av) >= 0) (here we restrict the 
function —21n Ap((av) | ^) to (av) > 0, whereas in previous MAGIC results Xp was also computed 
for negative (av) values, and the quoted limit was (av)svt = (av) 2 . 7 i — (av), whenever (av) < 0, 
which is a more conservative choice). 
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Figure 1: 95% CL upper limits on the thermally-averaged cross section for DM particles annihi¬ 
lating into bb (upper plots) and (lower plots). Thick solid lines show the combined limits 

obtained by merging the Fermi-LAT observations of 15 dSphs (left plots) or Segue 1 (right plots) 
with MAGIC observations of Segue 1. Dashed lines show the individual MAGIC (short dashes) 
and Fermi-LKF (long dashes) limits. J-factor statistical uncertainties are included. For the Segue 1 
results, we also show (thin-solid line) the combined limits assuming a fixed J-factor (no statistical 
uncertainties). The thin-dotted line, green and yellow bands show, respectively, the median and the 
symmetrical, two-sided 68% and 95% containment bands for the distribution of limits under the 
null {Hq : (av)=0) hypothesis (see main text for more details). The red-dashed-dotted line shows 
the thermal relic cross section from [24]. 

5. Discussion and Conclusions 

This work presents, for the first time, limits to the DM annihilation cross-section from a com¬ 
prehensive analysis of gamma-ray data of energies between 500 MeV and 10 TeV. Using a com¬ 
mon, homogeneous analysis approach (both in the applied statistical methods and in the determi¬ 
nation of the J-factors), we have combined the MAGIC observations of Segue I with Ferm/-LAT 
observations of 15 dSphs. This allowed the computation of meaningful global DM limits, and the 
direct comparison of the individual results obtained with different instruments. Our results span the 
DM particle mass range from 10 GeV to 100 TeV - the widest range covered by a single analysis 
to date. 

We have not observed any DM signal. Consequently, we set limits on the DM annihilation 
cross-section. Our results are the most constraining from observations of dSphs in the considered 
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mass range. For the low-mass range, our results (fully dominated by Fermi-LKT data) are below 
the thermal relie eross-seetion (av) ~ 3 x 10^^^ em^ s^'. In the intermediate mass range (from few 
hundred GeV to few tens TeV, depending on the eonsidered annihilation ehannel), where Fermi- 
LAT and MAGIC aehieve similar sensitivities, the improvement of the eombined result with respeet 
to the individual ones reaehes a faetor ~2. In addition, we present, for the first time, limits to high 
DM partiele mass above 10 TeV (fully dominated by MAGIC). 

Our global analysis method is eompletely generie, and ean be easily extended to inelude data 
from more targets, instruments and/or messenger partieles provided they have similar sensitivity to 
the eonsidered DM partiele mass range. Of partieular interest is the ease of a global DM seareh 
from dSphs ineluding data from all eurrent gamma-ray {Fermi-LPS, MAGIC, VERITAS, H.E.S.S, 
HAWC) and neutrino (leeCube, Antares) instruments, and we hereby propose a eoordinated effort 
toward that end. Ineluding results obtained from other types of observational targets like the Galae- 
tie Center, galaxy elusters or others is formally also possible, but a common approach to the J-factor 
determination remains an open question. In the future, this analysis could include new instruments 
like CTA, Gamma-400 or Km3Net. Our global approach offers the best chances for indirect DM 
discovery, or for setting the most stringent limits attainable by these kinds of observations, therefore 
placing a new landmark in the field. 
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